#'
#' quantify the dynamic change of phylogenetic or functional community structure
#' 
#' @param com1,com2 communities at different observation times
#' @param focus1i,focus2i the index of the same focus indiviudals in the first and second survey
#' @param distm phylogenetic or functional species pairwise distance matrix
#' @param r a vector of distances that the mark correlation change function should be evaluated
#' @param h a width of rectangle in the rectangle kernel smooth. see \link{define_rectSmooth} for more detials
#' 
#' @details
#' The mark correlation change function is defined as the mean neighborhood phylogenetic/functional structure changes
#' at scales between (r-h,r+h). it is designed to compare the dynamic of phylogenetic/functional structure during an
#' time interval. Thus positive values mean the structure becomes overdispersion, and negative values suggest the 
#' neighborhood structure becomes more clustering.
#' 
#' In order to test whether these none zero changes are significantly different with random community assembly,
#' a constrained random death and birth null model can be used (e.g. ... already implimented in this package). 
#' more sophisticated null models can be used here too to test more complex assumptions about community assembly.
#' 
#' Comparing with the phylogenetic mark correlation function, it does not need to be normalized by plot mean pairwise
#' phylogenetic/functional distance.
#' 
#' @return
#' a community structure object with mean and sd at various spatial scale r.
#' 
#' @examples
#' 
#' data(BCI)
#' BCI1=BCI
#' N=total_abundance(BCI)
#' focus1i=sample(1:N,1000)
#' alive=unique(c(focus1i,sample(1:N,20000)))
#' focus2i=match(focus1i,alive)
#' BCI2=subset(BCI,alive)
#' S=total_richness(BCI)
#' distm=as.matrix(dist(runif(S)))
#' r=seq(0,40,2)
#' h=1
#' re=MC_change(BCI1,BCI2,focus1i,focus2i,distm,r,h)
#' re
#' plot(re)

MC_change=function(com1,com2,focus1i,focus2i,distm,r,h){
  if(length(focus1i)!=length(focus2i)){
    stop("length of focus1i and focus2i should be the same")
  }else if( any(com1$species[focus1i]!=com2$species[focus2i]) ){
    stop("focus1i and focus2i not refer to the same individuals, \n Please reconsider the indexes of the same individuals in the given two community data")
  }
    
  rRange=range(r)
  dis1=frnn(com1,com1$x[focus1i],com1$y[focus1i],rRange=rRange,type="circle",info="both")
  dis2=frnn(com2,com2$x[focus2i],com2$y[focus2i],rRange=rRange,type="circle",info="both")
  
  Nf=length(focus1i)
  #change species name into index, and make sure the same index refers to the same species in community data and distm
  mcc=matrix(nrow=Nf,ncol=length(r))
  
  rdg=define_rectSmooth(r,h)
  #note that the factor level of species in com1 and com2 should be the same
  species_num1=as.numeric(com1$species)
  species_num2=as.numeric(com2$species)
  for(i in 1:Nf){
    sp1i=species_num1[focus1i[i]]
    sp2i=species_num2[focus2i[i]]
   ni1=meanDistr(species_num1,distm[sp1i,],dis1[[i]]$index,dis1[[i]]$dist,rdg) 
   ni2=meanDistr(species_num1,distm[sp2i,],dis2[[i]]$index,dis2[[i]]$dist,rdg) 
   mcc[i,]=ni1-ni2
  }
    
  result=data.frame(r=r,mean=apply(mcc,2,mean),sd=apply(mcc,2,sd))
  attr(result,"cname")="Mark correlation change function"
  class(result)=c("comstr",class(result))
  return(result)
}

meanDistr=function(species,distm,index,dist,rdg){
  
  dd=distm[species[index]]
  
  disti=apply_rectSmooth(rdg,dist)
  return(unlist(lapply(disti,function(x) ifelse(length(x)==0,0,mean(dd[x])) )))
}


plot.comstr=function(data,...){
  require(ggplot2)
  ggplot(data,aes(x=r,y=mean))+
    geom_line(aes(group=1))+
    geom_point(size=4)+
    geom_errorbar(aes(ymin=mean-sd,ymax=mean+sd))+
    xlab("Spatial distance")+
    ylab(attr(data,"cname"))
  
}
